home *** CD-ROM | disk | FTP | other *** search
/ Nebula 2 / Nebula Two.iso / Apps / Astro / ephem / Source / objx.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-06-12  |  28.8 KB  |  1,152 lines
  1. /* functions to save the user-definable objects, referred to as "x" and "y".
  2.  * this way, once defined, the objects can be quieried for position just like
  3.  * the other bodies, with obj_cir(). 
  4.  */
  5.  
  6. #include <stdio.h>
  7. #include <math.h>
  8. #include <ctype.h>
  9. #ifdef VMS
  10. #include <stdlib.h>
  11. #endif
  12. #ifdef NeXT
  13. # include <stdlib.h>
  14. #endif
  15. #include "astro.h"
  16. #include "circum.h"
  17. #include "screen.h"
  18.  
  19. extern char *strcat(), *strcpy(), *strncpy(), *getenv();
  20.  
  21. static char *dbfile;            /* !0 if set by -d option */
  22. static char dbfdef[] = "ephem.db";     /* default database file name */
  23.  
  24. /* structures to describe objects of various types.
  25.  */
  26. #define    MAXNM        16    /* longest allowed object name, inc \0 */
  27. typedef struct {
  28.     double m_m1, m_m2;    /* either g/k or H/G, depending on... */
  29.     int m_whichm;    /* one of MAG_gk or MAG_HG */
  30. } Mag;
  31. typedef struct {
  32.     double f_ra;    /* ra, rads, at given epoch */
  33.     double f_dec;    /* dec, rads, at given epoch */
  34.     double f_mag;    /* visual magnitude */
  35.     double f_siz;    /* angular size, in arc seconds */
  36.     double f_epoch;    /* the given epoch, as an mjd */
  37.     char   f_name[MAXNM]; /* name */
  38. } ObjF;            /* fixed object */
  39. typedef struct {
  40.     double e_inc;    /* inclination, degrees */
  41.     double e_Om;    /* longitude of ascending node, degrees */
  42.     double e_om;    /* argument of perihelion, degress */
  43.     double e_a;        /* mean distance, aka, semi-maj axis, in AU */
  44.     double e_n;        /* daily motion, degrees/day */
  45.     double e_e;        /* eccentricity */
  46.     double e_M;        /* mean anomaly, ie, degrees from perihelion at... */
  47.     double e_cepoch;    /* epoch date (M reference), as an mjd */
  48.     double e_epoch;    /* equinox year (inc/Om/om reference), as an mjd */
  49.     Mag    e_mag;    /* magnitude */
  50.     double e_siz;    /* angular size, in arc seconds at 1 AU */
  51.     char   e_name[MAXNM]; /* name */
  52. } ObjE;            /* object in heliocentric elliptical orbit */
  53. typedef struct {
  54.     double h_ep;    /* epoch of perihelion, as an mjd */
  55.     double h_inc;    /* inclination, degs */
  56.     double h_Om;    /* longitude of ascending node, degs */
  57.     double h_om;    /* argument of perihelion, degs. */
  58.     double h_e;        /* eccentricity */
  59.     double h_qp;    /* perihelion distance, AU */
  60.     double h_epoch;    /* equinox year (inc/Om/om reference), as an mjd */
  61.     double h_g, h_k;    /* magnitude model coefficients */
  62.     double h_siz;    /* angular size, in arc seconds at 1 AU */
  63.     char   h_name[MAXNM]; /* name */
  64. } ObjH;            /* object in heliocentric parabolic trajectory */
  65. typedef struct {
  66.     double p_ep;    /* epoch of perihelion, as an mjd */
  67.     double p_inc;    /* inclination, degs */
  68.     double p_qp;    /* perihelion distance, AU */
  69.     double p_om;    /* argument of perihelion, degs. */
  70.     double p_Om;    /* longitude of ascending node, degs */
  71.     double p_epoch;    /* reference epoch, as an mjd */
  72.     double p_g, p_k;    /* magnitude model coefficients */
  73.     double p_siz;    /* angular size, in arc seconds at 1 AU */
  74.     char   p_name[MAXNM]; /* name */
  75. } ObjP;            /* object in heliocentric parabolic trajectory */
  76.  
  77. typedef struct {
  78.     int  o_type;    /* current object type; see flags, below */
  79.     int  o_on;        /* !=0 while current object is active */
  80.     ObjF o_f;        /* the fixed object */
  81.     ObjE o_e;        /* the elliptical orbit object */
  82.     ObjH o_h;        /* the hyperbolic orbit object */
  83.     ObjP o_p;        /* the parabolic orbit object */
  84. } Obj;
  85.  
  86. /* o_type */
  87. #define    FIXED        1
  88. #define    ELLIPTICAL    2
  89. #define    HYPERBOLIC    3
  90. #define    PARABOLIC    4
  91.  
  92. /* m_whichm */
  93. #define    MAG_HG        0    /* using 0 makes HG the initial default */
  94. #define    MAG_gk        1
  95.  
  96. static Obj objx;
  97. static Obj objy;
  98.  
  99. #define    DY    0        /* decimal year flag for set_year() */
  100. #define    YMD    1        /* year/mon/day flag for set_year() */
  101.  
  102. /* run when Objx or y is picked from menu.
  103.  * we tell which by the planet code.
  104.  * let op define object and turn it on and off.
  105.  */
  106. obj_setup(p)
  107. int p;
  108. {
  109.     static char *pr[6] = { /* leave a slot for "On"/"Off" */
  110.         "Fixed", "Elliptical", "Hyperbolic", "Parabolic", "Lookup"
  111.     };
  112.     int f;
  113.     Obj *op;
  114.  
  115.     op = (p == OBJX) ? &objx : &objy;
  116.  
  117.     rechk:
  118.     /* map o_type to popup choice.
  119.      */
  120.     switch (op->o_type) {
  121.     case FIXED: f = 0; break;
  122.     case ELLIPTICAL: f = 1; break;
  123.     case HYPERBOLIC: f = 2; break;
  124.     case PARABOLIC: f = 3; break;
  125.     default: f = 4; break;
  126.     }
  127.  
  128.     ask:
  129.     pr[5] = op->o_on ? "On" : "Off";
  130.     switch (f = popup (pr, f, 6)) {
  131.     case 0: obj_dfixed(op, 0, (char**)0); goto ask;
  132.     case 1: obj_delliptical(op, 0, (char**)0); goto ask;
  133.     case 2: obj_dhyperbolic(op, 0, (char**)0); goto ask;
  134.     case 3: obj_dparabolic(op, 0, (char**)0); goto ask;
  135.     case 4: if (obj_filelookup(p, (char *)0) == 0) obj_on(p); goto rechk;
  136.     case 5: op->o_on ^= 1; break;
  137.     }
  138. }
  139.  
  140. /* turn "on" or "off" but don't forget facts about object the object.
  141.  */
  142. obj_on (p)
  143. int p;
  144. {
  145.     if (p == OBJX)
  146.         objx.o_on = 1;
  147.     else
  148.         objy.o_on = 1;
  149. }
  150. obj_off (p)
  151. int p;
  152. {
  153.     if (p == OBJX)
  154.         objx.o_on = 0;
  155.     else
  156.         objy.o_on = 0;
  157. }
  158.  
  159. /* return true if object is now on, else 0.
  160.  */
  161. obj_ison(p)
  162. int p;
  163. {
  164.     return ((p == OBJX) ? objx.o_on : objy.o_on);
  165. }
  166.  
  167. /* set an alternate database file name.
  168.  * N.B. we assume the storage pointed to by name is permanent.
  169.  */
  170. obj_setdbfilename (name)
  171. char *name;
  172. {
  173.     dbfile = name;
  174. }
  175.  
  176. /* fill in info about object x or y.
  177.  * most arguments and conditions are the same as for plans().
  178.  * only difference is that mag is already apparent, not absolute magnitude.
  179.  * this is called by body_cir() for object x and y just like plans() is called
  180.  * for the planets.
  181.  */
  182. obj_cir (jd, p, lpd0, psi0, rp0, rho0, lam, bet, siz, mag)
  183. double jd;    /* mjd now */
  184. int p;        /* OBJX or OBJY */
  185. double *lpd0;    /* heliocentric longitude, or NOHELIO  */
  186. double *psi0;    /* heliocentric latitude, or 0 if *lpd0 set to NOHELIO */
  187. double *rp0;    /* distance from the sun, or 0 */
  188. double *rho0;    /* true distance from the Earth, or 0 */
  189. double *lam;    /* apparent geocentric ecliptic longitude */
  190. double *bet;    /* apparent geocentric ecliptic latitude */
  191. double *siz;    /* angular size of object, arc seconds */
  192. double *mag;    /* APPARENT magnitude */
  193. {
  194.     Obj *op = (p == OBJX) ? &objx : &objy;
  195.  
  196.     switch (op->o_type) {
  197.     case FIXED: {
  198.         double xr, xd;
  199.         xr = op->o_f.f_ra;
  200.         xd = op->o_f.f_dec;
  201.         if (op->o_f.f_epoch != jd)
  202.         precess (op->o_f.f_epoch, jd, &xr, &xd);
  203.         eq_ecl (jd, xr, xd, bet, lam);
  204.  
  205.         *lpd0 = NOHELIO;
  206.         *psi0 = *rp0 = *rho0 = 0.0;
  207.         *mag = op->o_f.f_mag;
  208.         *siz = op->o_f.f_siz;
  209.         }
  210.         break;
  211.  
  212.     case ELLIPTICAL: {
  213.         /* this is basically the same code as pelement() and plans()
  214.          * combined and simplified for the special case of osculating
  215.          * (unperturbed) elements.
  216.          * inputs have been changed to match the Astronomical Almanac.
  217.          * we have added reduction of elements using reduce_elements().
  218.          */
  219.         double dt, lg, lsn, rsn;
  220.         double nu, ea;
  221.         double ma, rp, lo, slo, clo;
  222.         double inc, psi, spsi, cpsi;
  223.         double y, lpd, rpd, ll, rho, sll, cll;
  224.         double om;        /* arg of perihelion */
  225.         double Om;        /* long of ascending node. */
  226.         double e;
  227.         int pass;
  228.  
  229.         dt = 0;
  230.         sunpos (jd, &lsn, &rsn);
  231.         lg = lsn + PI;
  232.         e = op->o_e.e_e;
  233.  
  234.         for (pass = 0; pass < 2; pass++) {
  235.  
  236.         reduce_elements (op->o_e.e_epoch, jd-dt, degrad(op->o_e.e_inc),
  237.                 degrad (op->o_e.e_om), degrad (op->o_e.e_Om),
  238.                 &inc, &om, &Om);
  239.  
  240.         ma = degrad (op->o_e.e_M
  241.                 + (jd - op->o_e.e_cepoch - dt) * op->o_e.e_n);
  242.         anomaly (ma, e, &nu, &ea);
  243.         rp = op->o_e.e_a * (1-e*e) / (1+e*cos(nu));
  244.         lo = nu + om;
  245.         slo = sin(lo);
  246.         clo = cos(lo);
  247.         spsi = slo*sin(inc);
  248.         y = slo*cos(inc);
  249.         psi = asin(spsi);
  250.         lpd = atan(y/clo)+Om;
  251.         if (clo<0) lpd += PI;
  252.         range (&lpd, 2*PI);
  253.         cpsi = cos(psi);
  254.         rpd = rp*cpsi;
  255.         ll = lpd-lg;
  256.         rho = sqrt(rsn*rsn+rp*rp-2*rsn*rp*cpsi*cos(ll));
  257.         dt = rho*5.775518e-3;    /* light travel time, in days */
  258.         if (pass == 0) {
  259.             *lpd0 = lpd;
  260.             *psi0 = psi;
  261.             *rp0 = rp;
  262.             *rho0 = rho;
  263.         }
  264.         }
  265.  
  266.         sll = sin(ll);
  267.         cll = cos(ll);
  268.         if (rpd < rsn)
  269.         *lam = atan(-1*rpd*sll/(rsn-rpd*cll))+lg+PI;
  270.         else
  271.         *lam = atan(rsn*sll/(rpd-rsn*cll))+lpd;
  272.         range (lam, 2*PI);
  273.         *bet = atan(rpd*spsi*sin(*lam-lpd)/(cpsi*rsn*sll));
  274.  
  275.         if (op->o_e.e_mag.m_whichm == MAG_HG) {
  276.         /* the H and G parameters from the Astro. Almanac.
  277.          */
  278.         double psi_t, Psi_1, Psi_2, beta;
  279.         beta = acos((rp*rp + rho*rho - rsn*rsn)/ (2*rp*rho));
  280.         psi_t = exp(log(tan(beta/2.0))*0.63);
  281.         Psi_1 = exp(-3.33*psi_t);
  282.         psi_t = exp(log(tan(beta/2.0))*1.22);
  283.         Psi_2 = exp(-1.87*psi_t);
  284.         *mag = op->o_e.e_mag.m_m1 + 5.0*log10(rp*rho)
  285.             - 2.5*log10((1-op->o_e.e_mag.m_m2)*Psi_1
  286.             + op->o_e.e_mag.m_m2*Psi_2);
  287.         } else {
  288.         /* the g/k model of comets */
  289.         *mag = op->o_e.e_mag.m_m1 + 5*log10(rho)
  290.                     + 2.5*op->o_e.e_mag.m_m2*log10(rp);
  291.         }
  292.         *siz = op->o_e.e_siz / rho;
  293.         }
  294.         break;
  295.  
  296.     case HYPERBOLIC: {
  297.         double dt, lg, lsn, rsn;
  298.         double nu, ea;
  299.         double ma, rp, lo, slo, clo;
  300.         double inc, psi, spsi, cpsi;
  301.         double y, lpd, rpd, ll, rho, sll, cll;
  302.         double om;        /* arg of perihelion */
  303.         double Om;        /* long of ascending node. */
  304.         double e;
  305.         double a, n;    /* semi-major axis, mean daily motion */
  306.         int pass;
  307.  
  308.         dt = 0;
  309.         sunpos (jd, &lsn, &rsn);
  310.         lg = lsn + PI;
  311.         e = op->o_h.h_e;
  312.         a = op->o_h.h_qp/(e - 1.0);
  313.         n = .98563/sqrt(a*a*a);
  314.  
  315.         for (pass = 0; pass < 2; pass++) {
  316.  
  317.         reduce_elements (op->o_h.h_epoch, jd-dt, degrad(op->o_h.h_inc),
  318.                 degrad (op->o_h.h_om), degrad (op->o_h.h_Om),
  319.                 &inc, &om, &Om);
  320.  
  321.         ma = degrad ((jd - op->o_h.h_ep - dt) * n);
  322.         anomaly (ma, e, &nu, &ea);
  323.         rp = a * (e*e-1.0) / (1.0+e*cos(nu));
  324.         lo = nu + om;
  325.         slo = sin(lo);
  326.         clo = cos(lo);
  327.         spsi = slo*sin(inc);
  328.         y = slo*cos(inc);
  329.         psi = asin(spsi);
  330.         lpd = atan(y/clo)+Om;
  331.         if (clo<0) lpd += PI;
  332.         range (&lpd, 2*PI);
  333.         cpsi = cos(psi);
  334.         rpd = rp*cpsi;
  335.         ll = lpd-lg;
  336.         rho = sqrt(rsn*rsn+rp*rp-2*rsn*rp*cpsi*cos(ll));
  337.         dt = rho*5.775518e-3;    /* light travel time, in days */
  338.         if (pass == 0) {
  339.             *lpd0 = lpd;
  340.             *psi0 = psi;
  341.             *rp0 = rp;
  342.             *rho0 = rho;
  343.         }
  344.         }
  345.  
  346.         sll = sin(ll);
  347.         cll = cos(ll);
  348.         if (rpd < rsn)
  349.         *lam = atan(-1*rpd*sll/(rsn-rpd*cll))+lg+PI;
  350.         else
  351.         *lam = atan(rsn*sll/(rpd-rsn*cll))+lpd;
  352.         range (lam, 2*PI);
  353.         *bet = atan(rpd*spsi*sin(*lam-lpd)/(cpsi*rsn*sll));
  354.  
  355.         *mag = op->o_h.h_g + 5*log10(rho) + 2.5*op->o_h.h_k*log10(rp);
  356.         *siz = op->o_h.h_siz / rho;
  357.         }
  358.         break;
  359.  
  360.     case PARABOLIC: {
  361.         double inc, om, Om;
  362.         double lpd, psi, rp, rho;
  363.         double dt;
  364.         int pass;
  365.  
  366.         /* two passes to correct lam and bet for light travel time. */
  367.         dt = 0.0;
  368.         for (pass = 0; pass < 2; pass++) {
  369.         reduce_elements (op->o_p.p_epoch, jd-dt, degrad(op->o_p.p_inc),
  370.             degrad(op->o_p.p_om), degrad(op->o_p.p_Om), &inc, &om, &Om);
  371.         comet (jd-dt, op->o_p.p_ep, inc, om, op->o_p.p_qp, Om,
  372.                     &lpd, &psi, &rp, &rho, lam, bet);
  373.         if (pass == 0) {
  374.             *lpd0 = lpd;
  375.             *psi0 = psi;
  376.             *rp0 = rp;
  377.             *rho0 = rho;
  378.         }
  379.         dt = rho*5.775518e-3;    /* au to light-days */
  380.         }
  381.         *mag = op->o_p.p_g + 5*log10(rho) + 2.5*op->o_p.p_k*log10(rp);
  382.         *siz = op->o_p.p_siz / rho;
  383.         }
  384.         break;
  385.  
  386.     default:
  387.         f_msg ((p == OBJX) ? "Obj X is not defined"
  388.                    : "Obj Y is not defined");
  389.         break;
  390.     }
  391. }
  392.  
  393. /* define obj based on the ephem.db line, s.
  394.  * p is one of OBJX or OBJY.
  395.  * format: name,type,[other fields, as per corresponding ObjX typedef]
  396.  * N.B. we replace all ',' within s with '\0' IN PLACE.
  397.  * return 0 if ok, else print reason why not with f_msg() and return -1.
  398.  */
  399. obj_define (p, s)
  400. int p;    /* OBJX or OBJY */
  401. char *s;
  402. {
  403. #define    MAXARGS    20
  404.     char *av[MAXARGS];    /* point to each field for easy reference */
  405.     char c;
  406.     int ac;
  407.     Obj *op = (p == OBJX) ? &objx : &objy;
  408.  
  409.     /* parse into comma separated fields */
  410.     ac = 0;
  411.     av[0] = s;
  412.     do {
  413.         c = *s++;
  414.         if (c == ',' || c == '\0') {
  415.         s[-1] = '\0';
  416.         av[++ac] = s;
  417.         }
  418.     } while (c);
  419.  
  420.     if (ac < 2) {
  421.         char buf[NC];
  422.         if (ac > 0)
  423.         (void) sprintf (buf, "No type for Object %s", av[0]);
  424.         else
  425.         (void) sprintf (buf, "No fields in %s", s);
  426.         f_msg (buf);
  427.         return (-1);
  428.     }
  429.  
  430.     /* switch out on type of object - the second field */
  431.     switch (av[1][0]) {
  432.     case 'f':
  433.         if (ac != 6 && ac != 7) {
  434.         char buf[NC];
  435.         (void) sprintf(buf,
  436.             "Need ra,dec,mag,D[,siz] for fixed object %s", av[0]);
  437.         f_msg (buf);
  438.         return (-1);
  439.         }
  440.         obj_dfixed (op, ac, av);
  441.         break;
  442.  
  443.     case 'e':
  444.         if (ac != 13 && ac != 14) {
  445.         char buf[NC];
  446.         (void) sprintf (buf,
  447.             "Need i,O,o,a,n,e,M,E,D,H/g,G/k[,siz] for elliptical object %s",
  448.                                     av[0]);
  449.         f_msg (buf);
  450.         return (-1);
  451.         }
  452.         obj_delliptical (op, ac, av);
  453.         break;
  454.  
  455.     case 'h':
  456.         if (ac != 11 && ac != 12) {
  457.         char buf[NC];
  458.         (void) sprintf (buf,
  459.             "Need T,i,O,o,e,q,D,g,k[,siz] for hyperbolic object %s", av[0]);
  460.         f_msg (buf);
  461.         return (-1);
  462.         }
  463.         obj_dhyperbolic (op, ac, av);
  464.         break;
  465.  
  466.     case 'p':
  467.         if (ac != 10 && ac != 11) {
  468.         char buf[NC];
  469.         (void) sprintf (buf,
  470.             "Need T,i,o,q,O,D,g,k[,siz] for parabolic object %s", av[0]);
  471.         f_msg (buf);
  472.         return (-1);
  473.         }
  474.         obj_dparabolic (op, ac, av);
  475.         break;
  476.  
  477.     default: {
  478.         char buf[NC];
  479.         (void) sprintf (buf, "Unknown type for Object %s: %s",
  480.                                 av[0], av[1]);
  481.         f_msg (buf);
  482.         return (-1);
  483.         }
  484.     }
  485.  
  486.     return (0);
  487. }
  488.  
  489. /* if name, then look it up in the ephem database file and set p.
  490.  * else display a table of all objects and let op pick one.
  491.  * p is either OBJX or OBJY.
  492.  * if -d was used use it; else if EPHEMDB env set use it, else use default.
  493.  * return 0 if successfully set object p, else -1.
  494.  */
  495. obj_filelookup (p, name)
  496. int p;            /* OBJX or OBJY */
  497. char *name;
  498. {
  499. /* redefine RCTN,NTR,NTC for column-major order if you prefer */
  500. #define    NLR        (NR-1)        /* number of rows of names.
  501.                      * leave 1 for prompt
  502.                      */
  503. #define    LCW        9        /* screen columns per name */
  504. #define    NLC        9        /* total number of name columns */
  505. #define    NL        (NLR*NLC)    /* total number of names per screen */
  506. #define    RCTN(r,c)    ((r)*NLC+(c))    /* row/col to index */
  507. #define    NTR(n)        ((n)/NLC)    /* index to row */
  508. #define    NTC(n)        ((n)%NLC)    /* index to col (0 based) */
  509.                     /* N.B. all these are 0-based */
  510.     static char prompt[] =
  511.             "RETURN to select, p/n for previous/next page, q to quit";
  512.     FILE *fp;
  513.     char *fn;
  514.     int i, pgn;    /* index on current screen, current page number */
  515.     int r, c;
  516.     char buf[160];    /* longer than any one database line */
  517.     char pb[NC];    /* prompt buffer */
  518.     int readahd;    /* 1 if buffer set from previous loop */
  519.     int choice;    /* index to selection; -1 until set */
  520.     int roaming;    /* 1 while just roaming around screen */
  521.     int abandon;    /* 1 if decide to not pick afterall */
  522.  
  523.     /* open the database file */
  524.     if (dbfile)
  525.         fn = dbfile;
  526.     else {
  527.         fn = getenv ("EPHEMDB");
  528.         if (!fn)
  529.         fn = dbfdef;
  530.     }
  531.     fp = fopen (fn, "r");
  532.     if (!fp) {
  533.         (void) sprintf (buf, "Can not open database file %s", fn);
  534.         f_msg(buf);
  535.         return(-1);
  536.     }
  537.  
  538.     /* name is specified so just search for it without any op interaction */
  539.     if (name) {
  540.         int nl = strlen (name);
  541.         int ret = 0;
  542.         while (nxt_db(buf, sizeof(buf), fp) == 0 && strncmp(buf, name, nl))
  543.         continue;
  544.         if (feof(fp)) {
  545.         (void) sprintf (buf, "Object %s not found", name);
  546.         f_msg (buf);
  547.         ret = -1;
  548.         } else
  549.         (void) obj_define (p, buf);
  550.         (void) fclose (fp);
  551.         return (ret);
  552.     }
  553.         
  554.     pgn = 0;
  555.     readahd = 0;
  556.     choice = -1;
  557.     abandon = 0;
  558.  
  559.     /* continue until a choice is made or op abandons the attempt */
  560.     do {
  561.         /* put up next screen full of names.
  562.          * leave top row open for messages.
  563.          */
  564.         c_erase();
  565.         for (i = 0; i < NL; )
  566.         if (readahd || nxt_db (buf, sizeof(buf), fp) == 0) {
  567.             char objname[LCW];
  568.             int ii;
  569.             for (ii = 0; ii < sizeof(objname)-1; ii++)
  570.             if ((objname[ii] = buf[ii]) == ',')
  571.                 break;
  572.             objname[ii] = '\0';
  573.             if (i == NL-1)
  574.             objname[LCW-2] = '\0'; /* avoid scroll in low-r corner*/
  575.             f_string (NTR(i)+2, NTC(i)*LCW+1, objname);
  576.             i++;
  577.             readahd = 0;
  578.         } else
  579.             break;
  580.  
  581.         /* read another to check for eof. if valid, set readahd for next
  582.          * time.
  583.          */
  584.         if (nxt_db (buf, sizeof(buf), fp) == 0)
  585.         readahd = 1;
  586.  
  587.         /* let op pick one. set cursor on first one.
  588.          * remember these r/c are 0-based, but c_pos() is 1-based 
  589.          */
  590.         (void) sprintf (pb, "Page %d%s. %s", pgn+1,
  591.                         feof(fp) ? " (last)" : "", prompt);
  592.         f_prompt(pb);
  593.         r = c = 0;
  594.         roaming = 1;
  595.         do {
  596.         c_pos (r+2, c*LCW+1);
  597.         switch (read_char()) {
  598.         case 'h': /* left */
  599.             if (c == 0) c = NLC;
  600.             c -= 1;
  601.             if (RCTN(r,c) >= i)
  602.             c = NTC(i-1);
  603.             break;
  604.         case 'j': /* down */
  605.             if (++r == NLR) r = 0;
  606.             if (RCTN(r,c) >= i)
  607.             r = 0;
  608.             break;
  609.         case 'k': /* up */
  610.             if (r == 0) r = NLR;
  611.             r -= 1;
  612.             while (RCTN(r,c) >= i)
  613.             r -= 1;
  614.             break;
  615.         case 'l': /* right */
  616.             if (++c == NLC) c = 0;
  617.             if (RCTN(r,c) >= i)
  618.             c = 0;
  619.             break;
  620.         case REDRAW:
  621.             /* start over and skip over prior pages' entries */
  622.             rewind(fp);
  623.             for (i = 0; i < NL*pgn; i++)
  624.             (void) nxt_db (buf, sizeof(buf), fp);
  625.             readahd = 0;
  626.             roaming = 0;
  627.             break;
  628.         case 'p':
  629.             /* if not at first page, start over and skip back one
  630.              * pages' entries
  631.              */
  632.             if (pgn > 0) {
  633.             rewind(fp);
  634.             pgn--;
  635.             for (i = 0; i < NL*pgn; i++)
  636.                 (void) nxt_db (buf, sizeof(buf), fp);
  637.             readahd = 0;
  638.             roaming = 0;
  639.             }
  640.             break;
  641.         case 'n':
  642.             /* if not already at eof, we can go ahead another page */
  643.             if (!feof (fp)) {
  644.             pgn++;
  645.             roaming = 0;
  646.             }
  647.             break;
  648.         case END:
  649.             abandon = 1;
  650.             roaming = 0;
  651.             break;
  652.         case ' ': case '\r':
  653.             choice = NL*pgn + RCTN(r,c);
  654.             roaming = 0;
  655.             break;
  656.         }
  657.         } while (roaming);
  658.     } while (choice < 0 && !abandon);
  659.  
  660.     if (choice >= 0) {
  661.         /* skip first choice entries; selection is the next one */
  662.         (void) rewind (fp);
  663.         for (i = 0; i < choice; i++)
  664.         (void) nxt_db (buf, sizeof(buf), fp);
  665.         (void) nxt_db (buf, sizeof(buf), fp);
  666.         (void) obj_define (p, buf);
  667.     }
  668.     (void) fclose (fp);
  669.     redraw_screen (2);
  670.     return (choice >= 0 ? 0 : -1);
  671. }
  672.  
  673. /* read database file fp and put next valid entry into buf.
  674.  * return 0 if ok, else -1
  675.  */
  676. static
  677. nxt_db (buf, blen, fp)
  678. char buf[];
  679. int blen;
  680. FILE *fp;
  681. {
  682.     char s;
  683.     while (1) {
  684.         if (fgets (buf, blen, fp) == 0)
  685.         return (-1);
  686.         s = buf[0];
  687.         if (isalpha(s))
  688.         return (0);
  689.     }
  690. }
  691.  
  692. /* define a fixed object.
  693.  * args in av, in order, are name, type, ra, dec, magnitude, reference epoch
  694.  *   and optional angular size.
  695.  * if av then it is a list of strings to use for each parameter, else must
  696.  * ask for each (but type). the av option is for cracking the ephem.db line.
  697.  * if asking show current settings and leave unchanged if hit RETURN.
  698.  * END aborts without making any more changes.
  699.  * o_type is set to FIXED.
  700.  * N.B. we don't error check av in any way, not even for length.
  701.  */
  702. static
  703. obj_dfixed (op, ac, av)
  704. Obj *op;
  705. int ac;
  706. char *av[];
  707. {
  708.     char buf[NC];
  709.     char *bp;
  710.     int sts;
  711.  
  712.     op->o_type = FIXED;
  713.  
  714.     if (set_name (av, op->o_f.f_name) < 0)
  715.         return;
  716.  
  717.     if (av) {
  718.         bp = av[2];
  719.         sts = 1;
  720.     } else {
  721.         static char p[] = "RA (h:m:s): (";
  722.         f_prompt (p);
  723.         f_ra (R_PROMPT, C_PROMPT+sizeof(p)-1, op->o_f.f_ra);
  724.         (void) printf (") ");
  725.         sts = read_line (buf, 8+1);
  726.         if (sts < 0)
  727.         return;
  728.         bp = buf;
  729.     }
  730.     if (sts > 0) {
  731.         int h, m, s;
  732.         f_dec_sexsign (radhr(op->o_f.f_ra), &h, &m, &s);
  733.         f_sscansex (bp, &h, &m, &s);
  734.         sex_dec (h, m, s, &op->o_f.f_ra);
  735.         op->o_f.f_ra = hrrad(op->o_f.f_ra);
  736.     }
  737.  
  738.     if (av) {
  739.         bp = av[3];
  740.         sts = 1;
  741.     } else {
  742.         static char p[] = "Dec (d:m:s): (";
  743.         f_prompt (p);
  744.         f_gangle (R_PROMPT, C_PROMPT+sizeof(p)-1, op->o_f.f_dec);
  745.         (void) printf (") ");
  746.         sts = read_line (buf, 9+1);
  747.         if (sts < 0)
  748.         return;
  749.         bp = buf;
  750.     }
  751.     if (sts > 0) {
  752.         int dg, m, s;
  753.         f_dec_sexsign (raddeg(op->o_f.f_dec), &dg, &m, &s);
  754.         f_sscansex (bp, &dg, &m, &s);
  755.         sex_dec (dg, m, s, &op->o_f.f_dec);
  756.         op->o_f.f_dec = degrad(op->o_f.f_dec);
  757.     }
  758.  
  759.     if (set_double (av, 4, "Magnitude: ", &op->o_f.f_mag) < 0)
  760.         return;
  761.  
  762.     if (set_year (av, 5,"Reference epoch (UT Date, m/d.d/y or year.d): ",
  763.                             DY, &op->o_f.f_epoch) < 0)
  764.         return;
  765.  
  766.     if (ac == 7 || !av)
  767.         (void) set_double (av, 6, "Angular Size: ", &op->o_f.f_siz);
  768.     else
  769.         op->o_f.f_siz = 0.0;
  770.  
  771. }
  772.  
  773. /* define an object in an elliptical heliocentric orbit.
  774.  * 13 or 14 args in av, in order, are name, type, inclination, longitude of
  775.  *   ascending node, argument of perihelion, mean distance (aka semi-major
  776.  *   axis), daily motion, eccentricity, mean anomaly (ie, degrees from
  777.  *   perihelion), epoch date (ie, time of the mean anomaly value), equinox year
  778.  *   (ie, time of inc/lon/aop), two magnitude coefficients and optional size.
  779.  * the mag may be H/G or g/k model, set by leading g or H (use H/G if none).
  780.  * if av then it is a list of strings to use for each parameter, else must
  781.  * ask for each. the av option is for cracking the ephem.db line.
  782.  * if asking show current settings and leave unchanged if hit RETURN.
  783.  * END aborts without making any more changes.
  784.  * o_type is set to ELLIPTICAL.
  785.  * N.B. we don't error check av in any way, not even for length.
  786.  */
  787. static
  788. obj_delliptical(op, ac, av)
  789. Obj *op;
  790. int ac;
  791. char *av[];
  792. {
  793.     op->o_type = ELLIPTICAL;
  794.  
  795.     if (set_name (av, op->o_e.e_name) < 0)
  796.         return;
  797.  
  798.     if (set_double (av, 2, "Inclination (degs):", &op->o_e.e_inc) < 0)
  799.         return;
  800.  
  801.     if (set_double (av, 3, "Longitude of ascending node (degs):",
  802.                 &op->o_e.e_Om) < 0)
  803.         return;
  804.  
  805.     if (set_double (av, 4, "Argument of Perihelion (degs):",
  806.                 &op->o_e.e_om) < 0)
  807.         return;
  808.  
  809.     if (set_double (av, 5, "Mean distance (AU):", &op->o_e.e_a) < 0)
  810.         return;
  811.  
  812.     if (set_double (av, 6, "Daily motion (degs/day):", &op->o_e.e_n) < 0)
  813.         return;
  814.  
  815.     if (set_double (av, 7, "Eccentricity:", &op->o_e.e_e) < 0)
  816.         return;
  817.  
  818.     if (set_double (av, 8, "Mean anomaly (degs):", &op->o_e.e_M) < 0)
  819.         return;
  820.  
  821.     if (set_year (av, 9, "Epoch date (UT Date, m/d.d/y or year.d): ",
  822.                             YMD, &op->o_e.e_cepoch) < 0)
  823.         return;
  824.  
  825.     if (set_year (av, 10, "Equinox year (UT Date, m/d.d/y or year.d): ",
  826.                             DY, &op->o_e.e_epoch) < 0)
  827.         return;
  828.  
  829.     if (set_mag (av, 11, &op->o_e.e_mag) < 0)
  830.         return;
  831.  
  832.     if (ac == 14 || !av)
  833.         (void) set_double (av, 13, "Angular Size @ 1 AU: ", &op->o_e.e_siz);
  834.     else
  835.         op->o_e.e_siz = 0.0;
  836.  
  837. }
  838.  
  839. /* define an object in heliocentric hyperbolic orbit.
  840.  * 11 or 12 args in av, in order, are name, type, epoch of perihelion,
  841.  *   inclination, longitude of ascending node, argument of perihelion,
  842.  *    eccentricity, perihelion distance, reference epoch, absolute magnitude
  843.  *    and luminosity index, and optional size.
  844.  * if av then it is a list of strings to use for each parameter, else must
  845.  * ask for each. the av option is for cracking the ephem.db line.
  846.  * if asking show current settings and leave unchanged if hit RETURN.
  847.  * END aborts without making any more changes.
  848.  * o_type is set to HYPERBOLIC.
  849.  * N.B. we don't error check av in any way, not even for length.
  850.  */
  851. static
  852. obj_dhyperbolic (op, ac, av)
  853. Obj *op;
  854. int ac;
  855. char *av[];
  856. {
  857.     op->o_type = HYPERBOLIC;
  858.  
  859.     if (set_name (av, op->o_h.h_name) < 0)
  860.         return;
  861.  
  862.     if (set_year(av,2,"Epoch of perihelion (UT Date, m/d.d/y or year.d): ",
  863.                             YMD, &op->o_h.h_ep) < 0)
  864.         return;
  865.  
  866.     if (set_double (av, 3, "Inclination (degs):", &op->o_h.h_inc) < 0)
  867.         return;
  868.  
  869.     if (set_double (av, 4,
  870.         "Longitude of ascending node (degs):", &op->o_h.h_Om) < 0)
  871.         return;
  872.  
  873.     if (set_double(av,5,"Argument of perihelion (degs):", &op->o_h.h_om) <0)
  874.         return;
  875.  
  876.     if (set_double (av, 6, "Eccentricity:", &op->o_h.h_e) < 0)
  877.         return;
  878.  
  879.     if (set_double (av, 7, "Perihelion distance (AU):", &op->o_h.h_qp) < 0)
  880.         return;
  881.  
  882.     if (set_year (av, 8, "Reference epoch (UT Date, m/d.d/y or year.d): ",
  883.                             DY, &op->o_h.h_epoch) < 0)
  884.         return;
  885.  
  886.     if (set_double (av, 9, "g:", &op->o_h.h_g) < 0)
  887.         return;
  888.  
  889.     if (set_double (av, 10, "k:", &op->o_h.h_k) < 0)
  890.         return;
  891.  
  892.     if (ac == 12 || !av)
  893.         (void) set_double (av, 11, "Angular Size @ 1 AU: ", &op->o_h.h_siz);
  894.     else
  895.         op->o_h.h_siz = 0.0;
  896. }
  897.  
  898. /* define an object in heliocentric parabolic orbit.
  899.  * 10 or 11 args in av, in order, are name, type, epoch of perihelion,
  900.  *   inclination, argument of perihelion, perihelion distance, longitude of
  901.  *   ascending node, reference epoch, absolute magnitude and luminosity index,
  902.  *   and optional size.
  903.  * if av then it is a list of strings to use for each parameter, else must
  904.  * ask for each. the av option is for cracking the ephem.db line.
  905.  * if asking show current settings and leave unchanged if hit RETURN.
  906.  * END aborts without making any more changes.
  907.  * o_type is set to PARABOLIC.
  908.  * N.B. we don't error check av in any way, not even for length.
  909.  */
  910. static
  911. obj_dparabolic(op, ac, av)
  912. Obj *op;
  913. int ac;
  914. char *av[];
  915. {
  916.     op->o_type = PARABOLIC;
  917.  
  918.     if (set_name (av, op->o_p.p_name) < 0)
  919.         return;
  920.  
  921.     if (set_year(av,2,"Epoch of perihelion (UT Date, m/d.d/y or year.d): ",
  922.                             YMD, &op->o_p.p_ep) < 0)
  923.         return;
  924.  
  925.     if (set_double (av, 3, "Inclination (degs):", &op->o_p.p_inc) < 0)
  926.         return;
  927.  
  928.     if (set_double(av,4,"Argument of perihelion (degs):", &op->o_p.p_om) <0)
  929.         return;
  930.  
  931.     if (set_double (av, 5, "Perihelion distance (AU):", &op->o_p.p_qp) < 0)
  932.         return;
  933.  
  934.     if (set_double (av, 6,
  935.         "Longitude of ascending node (degs):", &op->o_p.p_Om) < 0)
  936.         return;
  937.  
  938.     if (set_year (av, 7, "Reference epoch (UT Date, m/d.d/y or year.d): ",
  939.                             DY, &op->o_p.p_epoch) < 0)
  940.         return;
  941.  
  942.     if (set_double (av, 8, "g:", &op->o_p.p_g) < 0)
  943.         return;
  944.  
  945.     if (set_double (av, 9, "k:", &op->o_p.p_k) < 0)
  946.         return;
  947.  
  948.     if (ac == 11 || !av)
  949.         (void) set_double (av, 10, "Angular Size @ 1 AU: ", &op->o_p.p_siz);
  950.     else
  951.         op->o_p.p_siz = 0.0;
  952. }
  953.  
  954. static
  955. set_double (av, vn, pr, fp)
  956. char *av[];    /* arg list */
  957. int vn;        /* which arg */
  958. char *pr;    /* prompt */
  959. double *fp;    /* ptr to double to be set */
  960. {
  961.     int sts;
  962.     char buf[NC];
  963.     char *bp;
  964.  
  965.     if (av) {
  966.         bp = av[vn];
  967.         sts = 1;
  968.     } else {
  969.         f_prompt (pr);
  970.         f_double (R_PROMPT, C_PROMPT+1+strlen(pr), "(%g) ", *fp);
  971.         sts = read_line (buf, 20);
  972.         if (sts < 0)
  973.         return (-1);
  974.         bp = buf;
  975.     }
  976.     if (sts > 0)
  977.         *fp = atof (bp);
  978.     return (0);
  979. }
  980.  
  981. static
  982. set_name (av, np)
  983. char *av[];    /* arg list */
  984. char *np;    /* name to be set */
  985. {
  986.     int sts;
  987.     char buf[NC];
  988.     char *bp;
  989.  
  990.     if (av) {
  991.         bp = av[0];
  992.         sts = 1;
  993.     } else {
  994.         (void) sprintf (buf, "Name: (%s) ", np);
  995.         f_prompt (buf);
  996.         sts = read_line (buf, MAXNM-1);
  997.         if (sts < 0)
  998.         return (-1);
  999.         bp = buf;
  1000.     }
  1001.     if (sts > 0)
  1002.         (void) strcpy (np, bp);
  1003.     return (0);
  1004. }
  1005.  
  1006. static
  1007. set_year (av, vn, pr, type, yp)
  1008. char *av[];    /* arg list */
  1009. int vn;        /* which arg */
  1010. char *pr;    /* prompt */
  1011. int type;    /* display type: YMD or DY */
  1012. double *yp;    /* ptr to year to be set */
  1013. {
  1014.     int sts;
  1015.     char buf[NC];
  1016.     char *bp;
  1017.  
  1018.     if (av) {
  1019.         bp = av[vn];
  1020.         sts = 1;
  1021.     } else {
  1022.         f_prompt (pr);
  1023.         if (type == DY) {
  1024.         double y;
  1025.         mjd_year (*yp, &y);
  1026.         (void) printf ("(%g) ", y);
  1027.         } else {
  1028.         int m, y;
  1029.         double d;
  1030.         mjd_cal (*yp, &m, &d, &y);
  1031.         (void) printf ("(%d/%g/%d) ", m, d, y);
  1032.         }
  1033.         sts = read_line (buf, 20);
  1034.         if (sts < 0)
  1035.         return (-1);
  1036.         bp = buf;
  1037.     }
  1038.     if (sts > 0)
  1039.         crack_year (bp, yp);
  1040.     return (0);
  1041. }
  1042.  
  1043. /* given either a decimal year (xxxx. something) or a calendar (x/x/x)
  1044.  * convert it to an mjd and store it at *p;
  1045.  */
  1046. static
  1047. crack_year (bp, p)
  1048. char *bp;
  1049. double *p;
  1050. {
  1051.     if (decimal_year(bp)) {
  1052.         double y = atof (bp);
  1053.         year_mjd (y, p);
  1054.     } else {
  1055.         int m, y;
  1056.         double d;
  1057.         mjd_cal (*p, &m, &d, &y);    /* init with current */
  1058.         f_sscandate (bp, &m, &d, &y);
  1059.         cal_mjd (m, d, y, p);
  1060.     }
  1061. }
  1062.  
  1063. /* read next two args from av and load the magnitude members m_m1 and m_m2.
  1064.  * also set m_whichm to default if this is from the .db file, ie, if av!=0.
  1065.  * #,#     -> model is unchanged
  1066.  * g#,[k]# -> g/k
  1067.  * H#,[G]# -> H/G
  1068.  */
  1069. static
  1070. set_mag (av, vn, mp)
  1071. char *av[];    /* arg list */
  1072. int vn;        /* which arg. we use av[vn] and av[vn+1] */
  1073. Mag *mp;
  1074. {
  1075.     int sts;
  1076.     char buf[NC];
  1077.     char *bp;
  1078.  
  1079.     if (av) {
  1080.         mp->m_whichm = MAG_HG;    /* always the default for the db file */
  1081.         bp = av[vn];
  1082.         sts = 1;
  1083.     } else {
  1084.         /* show both the value and the type of the first mag param,
  1085.          * as well as a hint as to how to set the type if desired.
  1086.          */
  1087.         (void) sprintf (buf, "%c: (%g) (g# H# or #) ",
  1088.                 mp->m_whichm == MAG_HG ? 'H' : 'g', mp->m_m1);
  1089.         f_prompt (buf);
  1090.         sts = read_line (buf, 9);
  1091.         if (sts < 0)
  1092.         return (-1);
  1093.         bp = buf;
  1094.     }
  1095.     if (sts > 0) {
  1096.         switch (bp[0]) {
  1097.         case 'g':
  1098.         mp->m_whichm = MAG_gk;
  1099.         bp++;
  1100.         break;
  1101.         case 'H':
  1102.         mp->m_whichm = MAG_HG;
  1103.         bp++;
  1104.         default:
  1105.         /* leave type unchanged if no prefix */
  1106.         break;
  1107.         }
  1108.         mp->m_m1 = atof (bp);
  1109.     }
  1110.  
  1111.     if (av) {
  1112.         bp = av[vn+1];
  1113.         sts = 1;
  1114.     } else {
  1115.         /* can't change the type in the second param */
  1116.         (void) sprintf (buf, "%c: (%g) ",
  1117.                 mp->m_whichm == MAG_HG ? 'G' : 'k', mp->m_m2);
  1118.         f_prompt (buf);
  1119.         sts = read_line (buf, 9);
  1120.         if (sts < 0)
  1121.         return (-1);
  1122.         bp = buf;
  1123.     }
  1124.     if (sts > 0) {
  1125.         int ok = 0;
  1126.         switch (bp[0]) {
  1127.         case 'k':
  1128.         if (mp->m_whichm == MAG_gk) {
  1129.             bp++;
  1130.             ok = 1;
  1131.         }
  1132.         break;
  1133.         case 'G':
  1134.         if (mp->m_whichm == MAG_HG) {
  1135.             bp++;
  1136.             ok = 1;
  1137.         }
  1138.         break;
  1139.         default:
  1140.         ok = 1;
  1141.         break;
  1142.         }
  1143.         if (ok)
  1144.         mp->m_m2 = atof (bp);
  1145.         else {
  1146.         f_msg ("Can't switch magnitude models at second parameter.");
  1147.         return (-1);
  1148.         }
  1149.     }
  1150.     return (0);
  1151. }
  1152.